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1. Introduction 

in 

CO Analysis in high energy physics and astrophysics is becoming increasingly 

complex. As part of the theoretical efforts tackling this complication, many 
multi-purpose event generators are developed. TASI program in 2011 pro- 
vided tutorials on selected tools during the school: CalcHEP, 1-4 PYTHIA, 5,6 
£SJ , MCFM, 7 PGS 8 and micrOmegas. 9,10 There are many other useful tools avail- 

able and users are strongly encouraged to experience those and compare 
with selected programs. Each tool has its advantages and disadvantages 
and users should be aware of their limitations. A Repository for Beyond- 
'*j_J ■ the-Standard- Model Tools can be found from several sources, some of which 

are listed below. 



(1) http:/ /www. ippp.dur.ac.uk/montecarlo/BSM/ 

(2) MC4BSM workshop: http://theory.fnal.gov/mc4bsm/ 

This note provides a quick summary for how to use CalcHEP and 
PYTHIA and emphasizes the use of batch modes that are often ig- 
nored. Advanced users should refer to manuals. 1-6 All plots and source 
codes are available from http://susy.phsx. ku.edu/^kckong/tasi/. Useful 
exercise for CalcHEP and PYTHIA are also found from PiTP school: 
http://www.phys. ufl.cdu/^matchev/pitp/ as well as their own web pages. 
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2. CalcHEP 

CalcHEP is a package for evaluation of Feynman diagrams, integration over 
multi-particle phase space, and event generation. It is a menu-driven system 
with contextual help. The notations used in CalcHEP are very similar to 
those used in particle physics. The CalcHEP package consists of two parts: 
symbolic and numerical. Both parts are written in the C programming 
language. The symbolic part produces C codes for a squared matrix element, 
and they are used in the numerical calculation later on. We summarize some 
features of CalcHEP below. 

• CalcHEP stands for Calculators for High Energy Physics. 

• CompHEP is a very similar program and shares many codes with 
CalcHEP 

— Download from http:/ /comphep. sinp.msu.ru/. 

— Current version is written in 2010. 

• Features 

— CalcHEP can evaluate any decay and scattering processes within 
any (user defined) model. 

— It has an easy user interface, and keeps correct spin correlations. 

— Symbolic calculation makes analysis very easy for 1 — > 2 and 2^2 
processes. 

— It is easy and quick to get plots from a model. 

— It is linked to mi cr Omegas for dark matter study. 

• Limitations 

— CalcHEP deals with tree-level processes only in the squared matrix 
element calculation. 

— It provides spin/polarization averaged amplitudes. 

— Limit on the number of external legs (involved particles) is 8 (for 
example, 2 — > 6 or 1 — > 7), and there is a limit on the number of 
diagrams. 

• CalcHEP is especially useful 

— when users need quick estimation of cross sections and decay 
widths (at tree- level), 

— when users cross-check calculations of other tools, 

— when users need to evaluate relic abundance of dark matter and 
direct /indirect detection limit, 

— when processes are relatively short. 
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2.1. Installation 

CalcHEP can be downloaded from 

http : //theory . sinp . msu . ru/~pukhov/calchep . html . 

To install, 

(1) unpack: tar -zxvf calchep_3.0.tar.gz 

(2) install by typing "make" . 

(3) If you encounter a problem with "blind mode", you need to install 
libxll-dcv. 

(4) create your own working directory: (e.g.) mkUsrDir WORK 

(5) to run, go to the WORK directory and simply execute "./calchep" 

The following subdirectories are created under the WORK directory: bin, 
models, results, tmp. Each model (under "models" directory) is defined 
in terms of 5 files. We will look into those in the next section. Other than 
main source codes, $CALCHEP/bin and $CALCHEP/utile contain useful 
scripts and codes. For more information see "README" in each directory. 
We will go over some useful routines with examples. 

CalcHEP is very easy to use. All selections can be made with "ENTER" 
or "RETURN" key. To go to the previous menu, simply press "ESC" key. 
If there are any questions, users can get real-time help by pressing "Fl" . 

2.2. A Closer look into model files: 

Particles, Vertices, Parameters, Constraints, Libraries 

A physics model is defined by a set of files with each containing 
a different aspect of the model: "prtcls#.mdl" for particle definition, 
"vars#.mdl" for independent parameters, "func#.mdl" for dependent pa- 
rameters, "lgrng#.mdl" for interactions, and '"extlib#.mdl" for external 
routines. 

(1) "Particles" : particles arc defined in this file. For example, particle def- 
inition for the Standard Model is shown below (usually "prtclsl.mdl"). 

Standard Model 
Particles 

Full name I >A < I >A < I number I 2*spin I mass I width I color I aux I >LaTex(A)< I >LaTeX(A+) <l 



gluon 


IG 


IG 


121 


12 


10 


10 


18 


IG 


lg 


lg 


photon 


IA 


IA 


122 


12 


10 


10 


11 


IG 


1 \gamma 


1 \gamma 


Z-boson 


IZ 


IZ 


123 


12 


IMZ 


IwZ 


11 


IG 


IZ 


IZ 


W-boson 


IW+ 


IV- 


124 


12 


IMW 


IwW 


11 


IG 


|V+ 


IV- 


Higgs 


h 


Ih 


125 


10 


IMh 


1 !wh 


11 


1 


In 


Ih 


electron 


le 


IE 


111 


11 


IMe 


10 


11 


1 


le~- 


le"+ 


e-neutrino 


1 ne 


Ne 


112 


11 


10 


10 


11 


|L 


1 \nu_e 


1 \bar{\nu>_e 
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muon 


Im 


IM 


113 


11 


|Mm 


10 


11 


1 l\mu~- 


1 \mu~+ 


m-neutrino 


1 nm 


|Nm 


114 


11 


10 


10 


11 


1 L 1 \nu_\mu 


1 \bar{\nu}_\mu 


tau-lepton 


11 


IL 


115 


11 


1 Ml 


10 


11 


1 l\tau~- 


1 \tau~- 


t-neutrino 


Inl 


INI 


116 


11 


10 


10 


11 


IL l\nu_\tau 


1 \bar{\nu}_\tau 


d- quark 


Id 


ID 


11 


11 


10 


10 


13 


1 Id 


1 \bar{d} 


u- quark 


In 


IU 


12 


11 


10 


10 


13 


1 In 


1 \bar{u} 


s -quark 


Is 


IS 


13 


11 


10 


10 


13 


1 Is 


1 \bar{s} 


c -quark 


Ic 


IC 


14 


11 


IMc 


10 


13 


1 Ic 


1 \bar{c} 


b- quark 


lb 


IB 


15 


11 


1Mb 


10 


13 


1 lb 


1 \bar{b} 


t -quark 


It 


IT 


16 


11 


IMt 


Iwt 


13 


1 It 


1 \bar{t} 



Each column is self-explanatory, and has information about particle 
name, symbol for particle, symbol for its anti-particle, PDG num- 
ber, spin, mass, width, color charge, auxiliary field and latex expres- 
sion. All variables need to be defined either in "Parameters" or in 
"Constraints" except for a width that is prefixed with an exclama- 
tion mark "!" . For instance, the Higgs width is defined as "!wh" , which 
means it is calculated on the fly whenever needed. CalcHEP supplies 
its value by running a separate decay process and the variable does 
not need to be defined in other model files such as "Parameters" or 
"Constraints". 

(2) "Parameters" contains all necessary independent parameters that de- 
fine a model and a numerical value is assigned to each parameter. 
Users can change values of those parameters, when needed. Below is 
"varsl.mdl" for the Standard Model. 

Parameters 



>Name 


< 1 Value 


1 > Comment 


<l 


alfEMZ 


10. 0078180608 IMS-BAR electron^ 


;netic alpha(MZ) 


alfSMZ 


10.1172 


ISrtong alpha(MZ) 


for running mass calculation 


GG 


11.238 


1 Running Strong coupling. The given value doesn't matter. 


SW 


10.481 


IMS-BAR sine of the electroweak mixing angle 


Mm 


10.1057 


Imuon mass 




Ml 


11.777 


1 tau-lepton mass 




McMc 


11.2 


IMc(Mc) 




MbMb 


14.25 


1Mb (Mb) 




Mtp 


1 175 


It-quark pole mas 


5 


MZ 


191.187 


1 Z-boson mass 




Mh 


1 120 


1 higgs mass 




wt 


11.59 


It-quark width 


(tree level l->2x) 


wZ 


12.49444 


1 Z-boson width 


(tree level l->2x) 


WW 


12.08895 


IW-boson width 


(tree level l->2x) 



All other variables are calculated based on above parameters. 

(3) "Constraints" defines variables which depend on parameters or con- 
straints that are defined previously. Their values are calculated auto- 
matically when CalcHEP is running, "funcl.mdl" is shown below. 

Standard Model 
Constraints 
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>Name 


<|> Expression 




EE 


1 sqrt (16*atan(l . ) *alf EMZ) 


% electromagnetic constant 


CW 


Isqrt(l-SIT 2) 


7, cos of the Weinberg angle 


MW 


IMZ*CW 


% W-boson mass 


Cl2 


Isqrt(l-sl2~ 2) 


% parameter of C-K-M matrix 


c23 


Isqrt(l-s23~ 2) 


% parameter of C-K-M matrix 


cl3 


Isqrt(l-sl3~ 2) 


% parameter of C-K-M matrix 


Vud 


Icl2*cl3 


X C-K-M matrix element 



qcdOk I initQCD (alf SMZ ,McMc , MbMb , Mtp) 

Mb IMbEff (Q)*one(qcd0k) 

Mt IMtEff (Q)*one(qcd0k) 

Mc IMcEff (Q)*one(qcd0k) 



Basic arithmetic operations such as +, -, /, *, , SqrtQ are allowed. 
Also external C functions are also allowed and they need to be defined 
in "usrfun.c". 

(4) "Vertices" contains actual interactions of particles ("lgrngl.mdl"). 
Current version of CalcHEP allows interactions with 4 particles. An 
interaction is defined by a prefactor and Lorentz structure as shown 
below. Interestingly, in the Lorentz part, "/" is not allowed and so 
any division should be carried out in the "Factor" column. The m2.m3 
(m3.m4) represents the metric g m2 m 3 , where m2 (m3) is the Lorentz 
index of the second (third) particle. 



Standard Model 
Vertices 

Al IA2 I A3 IA4 l> Factor <l> Lorentz part 



h 


IW+ 


|W- 


1 


1 EE*MW/SW 


1 m2 . m3 


h 


IZ 


|Z 


1 


IEE/(SW*CW~ 2)*MW 


1 m2 . m3 


h 


In 


h 


1 


|-(3/2)*EE*Mh~ 2/(MW*SW) 


11 


h 


In 


In 


In 


1 (-3/4)*(EE*Mh/(MW*SW))~ 2 


11 


h 


In 


IZ 


IZ 


I (1/2)*(EE/(SW*CW))~ 2 


1 m3 . m4 


h 


In 


IW+ 


IW- 


1 (1/2)*(EE/SW)" 2 


1 m3 . m4 


M 


Im 


In 


1 


|-EE*Mm/(2*MW*SW) 


11 


L 


11 


In 


1 


|-EE*M1 /(2*MW*SW) 


11 


C 


Ic 


In 


1 


|-EE*Mc/(2*MW*SW) 


11 


B 


lb 


h 


1 


|-EE*Mb/(2*MW*SW) 


11 


T 


It 


In 


1 


1 -EE*Mt /(2*MW*SW) 


11 



The usual i that appears in the Feynman rules is omitted by default. 

(5) "Libraries" contains user-defined codes ("extlibl.mdl" is shown as an 
example below.). 

Standard Model 
Libraries 

External libraries and citation <l 
$CALCHEP/utile/usrfun.c 
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2.3. Running examples 

2.3.1. Tree level branching fractions of the Higgs 

To warm up we start with an easiest example: two body decay of the Higgs 
in the Standard Model. Now run CalcHEP and follow the steps below to 
open a numerical session. 

(1) Choose "Standard Model" in the main menu and turn on unitary gauge. 

(2) Select "Enter Process". CalcHEP shows the previous process on screen. 

(3) Input the two body decay of the Higgs as follows. We will not exclude 
any particles or any diagrams. So leave them as blanks. 

Enter process: h -> 2*x 
Exclude diagrams with 
Exclude X-particles 

(4) CalcHEP checks whether the directory "results/" is empty. Users can 
cither delete or rename results from the previous run. 

(5) In the next screen users can verify diagrams by choosing "view dia- 
grams" . Users can choose certain diagrams, if necessary. 

(6) To continue to a numerical session, select "Squaring technique" and 
then "Make & launch n_calchep" 

(7) New GUI window will appear with different menus. 

Once numerical session is opened, users can change parameters and 
impose cuts. For detailed information regarding each menu, users should 
refer to the manual. In this exercise, we simply use the last menu "Easy 
1— >2". CalcHEP should show the total width with branching fractions on 
the next screen. For a different Higgs mass (or other parameters in general) , 
users can play with "parameter dependence" in the menu. CalcHEP allows 
users to plot quantities, branching fractions and total width, for instance, 
as a function of one parameter in the model, and output the data into a 
file. 

To get the total decay width, we can also run a script from the command 
line. Quit current numerical session by pressing Esc key a few times. To use 
scripts, go to the "results" directory and run "subproc_cycle" as follows. 

cd results 

. ./bin/subproc_cycle 
width (h) =0 . 00303927735 

#Subprocess 2 ( h -> b, B ) width=0 . 0025155 Br=0 . 8276638524 
#Subprocess 3 ( h -> c, C ) width=0 . 0001038 Br=0 . 03415285545 
#Subprocess 4 ( h -> 1, L ) width=0 . 0002501 Br=0 . 08228929815 
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#Subprocess 5 ( h -> m, M ) width=8.8605E-07 Br=0 . 0002915331172 
#Subprocess 8 ( h -> A, Z ) width=9 . 049E-07 Br=0 . 0002977352495 
#Subprocess 9 ( h -> A, A ) width=5 . 9964E-06 Br=0 . 001972969002 
#Subprocess 10 ( h -> G, G ) width=0 . 00016209 Br=0 . 05333175664 

Here the input "0" after "subproc_cycle" is the total number of events. 
Output includes the total decay width as well as partial widths and branch- 
ing fractions a . The numerical session can be called again anytime by run- 
ning n_calchep which is in "results" directory. 

2.3.2. 4 body decay: h — > e~ + v e + ji + + 

Now let us look at a slightly more complicated exercise with the four body 
decay of the Higgs to two charged leptons via two Ws. We follow similar 
steps as before, but now input a different process 

h -> e, Ne, M, nm 

There should be two diagrams, one with two Ws and the other with one 
W. The latter appears due to the Yukawa coupling of the muon, which is 
negligible. Launch numerical session. Now the last menu in the previous 
run, "Easy 1 — >2", does not appear in this example. To calculate the four 
body partial decay width, users need to run "Vegas", where the number 
of iterations and the number of MC points are set up. To calculate the 
width, simple choose "Start integration" and the result will be shown on 
screen. To obtain a partial width for a different Higgs mass, users should 
provide a different numerical value to the mass variable in the menu "Model 
parameters" and run Vegas again. 

Alternative way to compute it is to use a batch mode b . To use it, quit 
numerical session first and go to result directory and run "name.cycle" that 
is in the "bin" directory. 

. . /bin/name_cycle Mh 50 10 10 

This script computes the corresponding decay for 10 different Higgs masses 
starting from 50 GeV with 10 GeV interval. Output is shown on screen and 
saved into a file as well. 

Mh=50 sigma=3.2620E-10[pb] 



a The last two decay modes are added through an exercise in this note. If you haven't 
done so, you may not see them in the output. 

b For the most recent version of the batch mode, see the manual. 4 
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Mh=60 sigma=1.3323E-09[pb] 
Mh=70 sigma=4.6880E-09[pb] 
Mh=80 sigma=1.4999E-08[pb] 
Mh=90 sigma=5.6488E-08[pb] 
Mh=100 sigma=3.4003E-07[pb] 
Mh=110 sigma=1.6884E-06[pb] 
Mh=120 sigma=6.2347E-06[pb] 
Mh=130 sigma=1.8358E-05[pb] 
Mh=140 sigma=5.0125E-05[pb] 
See Mh_tab_2_ll file. 



Taking the data file with more points, one can plot the four body decay 
width of the Higgs in the e~fi + final state, which is shown in Fig. 1 c . One 
can notice the slight change in the slopes near Mw and 2Myy- The former 
is due to the fact that one of the W becomes onshell, while in the latter 
both W are onshell. 



c Users are supposed to use their own plotting package for this. 
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2.3.3. pp -+ Wbb 

In this example, we consider the scattering process pp — > Wbb at the Teva- 
tron. Steps are given as follows. 

(1) Open numerical session with the following process. We will call a proton 
as "p" and will include the first generation quarks and gluon in it. 

Enter process: p, p -> W+, b, B 
composite 'p' consists of: u, d, U, D, G 
Exclude diagrams with 

(2) In "IN State", choose PDF sets and set up the momentum of each 
beam. For instance, CTEQ5M(proton) for "S.F.I" and CTEQ5M(anti- 
proton) "S.F.2", and 980 GcV for each beam. 

(3) Cuts, kinematics and regularizations are useful for fast and more accu- 
rate computation. An example for cuts is shown below. 

#Cuts 
*** Table *** 
Cuts 

! I Parameter I > Min bound < I > Max bound < I 



There are many other available variables that are used as cuts (get help 
by pressing Fl.): A (angle in degree units), C (cosine of angle), J (jet 
cone angle, which is defined as y/ Ay 2 + Ar/> 2 , where Ay is the pseudo- 
rapidity difference and A<f> is the azimuthal angle difference between 
two particles), E (energy of the particle set), M (mass of the particle 
set), P (cosine in the rest frame of pair), T (transverse momentum of 
the particle set), S (squared invariant mass of the particle set), Y (ra- 
pidity of the particle set), and U (user's implemented function). The 
character string following U is passed on to the user C function usr- 
fun(str) which, as it is assumed, calculates a corresponding value. See 
Section 2.3 for further explanations. 



1Kb) 
IT(B) 
lN(b) 
IN(B) 
I J(b,B) 



120 
120 
1-5 
1-5 
10.4 



15 
15 



To set up relevant kinematics and regularizations, users need to take 
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a look at diagrams and specify them. For instance, in this example, b 
and b may be combined. 

#K inemat i c al _ s cheme 
12 -> 45 , 3 
45 -> 4 , 5 

The corresponding numbers are from the name of the process in the top 
line of numerical session of CalcHEP . They are numbered consecutively. 
As one can see from the process label, 

(sub)Process : u, D -> W+, b, B 

CalcHEP assigns 3 for W+, 4 for b and 5 for B. 

In the regularizations users specify a list of dangerous propagators. 
CalcHEP performs regularizations of the squared matrix element ac- 
cording to the contents of this list. 



#Regularization 
*** Table *** 
Regular izat ion 



Momentum 


1 > Mass 


<|> Width < 


1 


45 


IMZ 


1 wZ 


12 


45 


1Mb. 


Iwh 


12 


34 


iMt 


Iwt 


12 


35 


|Mt 


Iwt 


12 



To see distributions in CalcHEP , they need to be set up before VEGAS 
is called. Examples are list below. 

Parameter_l |> Min_l <|> Max_l <| 

M(W+, b) 80 200 

M(W+, B) 80 200 

M(b,B) 200 



Once users set up distributions, click on "Start integration" to get the 
total cross sections. For plots, select "Display Distributions". To im- 
prove errors and chi**2, increase nSess_l and nCalls_l (and nSess_2 
and nCalls_2 for 2 dimensional plots). 
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2.3.4. tt production cross section at the Tevatron 
Steps for pp — > tt are summarized as follows. 

(1) Run pp — > tt similarly with minor change as follows. 
Enter process: p, p -> t, T 

composite 'p' consists of: u, d, s, U, D, S, G 
Exclude diagrams with A, W+, W-, Z, h 

The last step forces CalcHEP to include diagrams which are mediated 
by strong interaction only. 

(2) Open numerical session and turn on PDF with CTEQ6L (for proton 
and anti-proton) and set the momentum of beams. This time we set 
QCD scale to Mt/2 (rrii 0p /2), which is known to minimize the difference 
between leading order and next-leading order cross sections. All this 
information is saved in files, "prt_#" and "session.dat" . 

#Initial_state inPl=9 . 800000E+02 inP2=9 . 800000E+02 
Polarizations= { 0.000000E+00 0.000000E+00 } 
StrFunl="PDT: cteq61 (proton) " 2212 
StrFun2="PDT: cteq61 (anti-proton) " -2212 

and 



#QCD alphaPDF=l alpha(MZ)=l . 172000E-01 NF=5 0rder=2 
MbMb=4.200000E+00 Mtp=l . 750000E+02 Scale= Mt/2 

(3) Use "subproc_cycle" to get the total cross section. 



. /bin/subproc_cycle 10 



#Subprocess 1 ( u, 

#Subprocess 2 ( d, 

#Subprocess 3 ( s, 

#Subprocess 4 ( U, 

#Subprocess 5 ( D, 

#Subprocess 6 ( S, 

#Subprocess 7 ( G, 



Cross section = 
Cross section = 
Cross section = 
Cross section = 
Cross section = 
Cross section = 
Cross section = 
1. 



Sum of distributions is stored in file distr 
Total Cross Section 7.8389739 [pb] 
see details in prt_l - prt_7 files 



1934E+00, 
1791E+00, 
1461E-03, 
3220E-02, 
9288E-02, 
1498E-03, 
2B67E-01, 



61934 events 
11791 events 
41 events 
132 events 
192 events 
41 events 
4256 events 



Total cross section summed over all subprocesses is reported at the end 
of output. The parameter that follows "subproc.cycle" is mandatory 
and is the total integrated luminosity in [1/fb] unit. Events are not 
generated with "subproc_cycle 0" . 
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(4) Now we can calculate the top mass dependence in the total cross sec- 
tions using "name_cycle_plus" d . 

. . /bin/name_cycle_plus Mtp 165 1 11 
Mtp=165 sigma=9. 6448496 [pb] 
Mtp=166 sigma=9. 3545769 [pb] 
Mtp=167 sigma=9. 0866098 [pb] 
Mtp=168 sigma=8. 8071721 [pb] 
Mtp=169 sigma=8. 5578445 [pb] 
Mtp=170 sigma=8. 3177445 [pb] 
Mtp=171 sigma=8. 064445 [pb] 
Mtp=172 sigma=7. 843451 [pb] 
Mtp=173 sigma=7. 6188025 [pb] 
Mtp=174 sigma=7. 3939838 [pb] 
Mtp=175 sigma=7. 1815069 [pb] 
See Mtp_cycle file. 

As expected from above output, "name_cycle_plus" reports the produc- 
tion cross section for 11 different masses with 1 GeV interval, starting 
from 165 GeV. 

(5) If distributions are defined before "subproc_cycle" is called, CalcHEP 
generates the corresponding distributions for each subprocess. Users 
can see them simply typing 

show_distr distr_# 

To combine distributions for all subprocess, run 
sum_distr distr_l distr_2 > distr_sum 

and run "show_distr distr_sum" for view. 

(6) "event_mixer" combines generated events files and write output in LHE 
format, that can be then used in any other event generators. It needs 
two inputs: the number of events and a directory where events files 
exist. If current path is 'results' directory, simply type 

. . /bin/event_mixer 10000 . 
total cross section 7.175E+00 
Max number of events 71503 



d This script is written by A. Pukhov but is not included in the CalcHEP package. Users 
need to download from http://susy.phsx.ku.edu/~kckong/tasi/name_cycle_plus. 
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and this combines all events in 7 subprocess and weight them ac- 
cording to their individual cross sections. The default output file is 
"event_mixer.lhe" and users can generate simple distributions with this 
LHE file, using "lhc2tab" . For example, try the following exercise. 

. ./bin/lhe2tab "M(6,-6)" 300 1000 100 < event_mixer . lhe > Mtt.txt 

The generated "Mtt.txt" files contains differential cross section as a 
function of the invariant mass of top (6) and anti-top (-6). 

2.3.5. User defined variables, and constraints in external C file 

User programs in CalcHEP provide users with possibility to attach his/her 
own codes to the n_calchep. In this way users are able to expand the set of 
phase space functions for cuts and histograms and implement new structure 
functions. We will discuss three examples: transverse mass, Mti and the 
Higgs effective coupling. We define the first two quantities in "usrfun.c" 
that is found under "utile" directory. An example is included as shown 
below. 

double usrfun(char * name, double * pvect) 
{ 

if (strcmp(name, "MT")==0) { /* Transverse mass */ 

double tmass, ET1 , ET2; 
double ppl [4] ; 
double pp2 [4] ; 
int k; 

for(k=0; k<4; k++) 
{ 

ppl[k] = pvect [4* (3-1) +k] ; 
pp2[k] = pvect [4* (4-1) +k] ; 

} 

ET1 = sqrt( pow(ppl [1] ,2) + pow(ppl [2] ,2) ); 
ET2 = sqrt( pow(pp2 [1] ,2) + pow(pp2 [2] ,2) ); 

tmass = 2 * ( ET1 * ET2 - ppl[l] * pp2[l] - ppl [2] * pp2 [2] ); 
return sqrt (tmass); 

} 

if (strcmp(name,"MT2")==0) { /* MT2 for massless particles */ 
double tmass2,ETl, ET2; 
double ppl [4] ; 
double pp2 [4] ; 
int k; 
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for(k=0; k<4; k++) 
{ 

ppl[k] = pvect[4*(3-l)+k] ; 
pp2[k] = pvect[4*(5-l)+k] ; 

} 



ET1 = sqrt( pow(ppl [1] , 2) + pow(ppl [2] ,2) ); 
ET2 = sqrt( pow(pp2 [1] , 2) + pow(pp2 [2] , 2) ); 

tmass2 = 2 * ( ET1 * ET2 + ppl [1] * pp2 [1] + ppl [2] * pp2 [2] ) ; 
return sqrt (tmass2) ; 

> 

/* original USRFUNC */ 

fprintf (stdout, " usrfun(char* name) called with parameter °/,s\n" 

" But is not def ined! \n" .name) ; 

sortie (54) ; 
return . ; 



Users also need to specify its location in a model file, "extlib9.mdl" . 

Standard Model 
Libraries 



$CALCHEP/utile/usrf un . c 



Now let us consider single production of the W and its leptonic decay. 
A well known and useful quantity in this process is the transverse mass, 
which is not implemented in CalcHEP . 

(1) Run p, p — >e, Ne at the Tevatron (or at the LHC). 

(2) To take a look at the invariant mass and the transverse mass of e and 
v e , we defined distributions as follows. 

Distributions 

Parameter_l I > Min_l <|> Max_l <| 



} 



External libraries and citation 



<l 



M(e,Ne) 
UMT 



I 
I 



100 I 
100 I 
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Here new distribution "UMT" is one of the functions that we defined 
in "usrfun.c" above. For CalcHEP to understand that this is a user 
defined function, the letter "U" should come with the name of user 
variables. This can be used in "cuts" as well. Unfortunately CalcHEP 
currently does not support having the particle names inside user-defined 
functions. 

Another good example of user-defined routines is the Higgs effective 
couplings to two photons. It is well known that the dominant production of 
the Higgs at the LHC is via glue-glue fusion. The main contribution arises at 
loops and unfortunately many tree- level event generators including CalcHEP 
do not treat this interaction. 

However one can parameterize them in term of higher dimensional op- 
erators and often the coefficients of such operators involve complicated ex- 
pressions such as integrals. In this case, one can perform calculation in a 
separate code, "usrfun.c" and return the result into CalcHEP . We will come 
back to this later in section 2.4.1. For discussion in the rest of this section, 
suppose we already have this interaction implemented. 

Next example is single production of the Higgs and its decay to two 
leptons via two Ws. 

(1) To Run gg -> ft ->• W + W~ -> e _ P e M +I/ A" type 

Enter process: G, G -> e, Ne , M, nm 
Exclude diagrams with m 

(2) Kinematics, regularizations and distributions are set up as follows. 

#Kinematical_scheme 
12 -> 34 , 56 
34 -> 3 , 4 
56 -> 5 , 6 

#Regularization 
Regular izat ion 



Momentum 


1 > Mass 


<|> Width 


<| Power I 


34 


IMW 


1 wW 


12 | 


56 


IMW 


1 wW 


12 | 


3456 


|Mh 


Iwh 


12 | 



#Distributions 
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Parameter_l 



l> 



Min_l <|> 
0.000000E+00I 
0.000000E+00I 
0.000000E+00I 
0.000000E+00I 



Max_l <| 
1.000000E+02I 
1.000000E+02I 
2.000000E+02I 
1.000000E+02I 



M(e,Ne) 
M(M,nm) 
M12 



UMT2 



The "M12" distribution is the invariant mass of the first two particles, 
i.e., G and G. Therefore this is equivalent to and the invariant mass 
of all particles in the final state, which should peak at the Higgs mass 
in this example. The "UMT2" is the Mti distribution taking e~ and 
\i + as visible particles and the two neutrinos as the missing particles. 
Once they are set up, perform Vegas integration, which will report cross 
sections at the end of the run. To see previously-defined distributions, 
click on "Display Distributions" . 

2.3.6. KK photon annihilation in mathematica 

One of excellent features of CalcHEP is that it provides analytic expres- 
sions for squared matrix elements in FORM, REDUCE and Mathematica. 
To learn how to exploit this analytic feature, we consider KK photon an- 
nihilation in Universal Extra Dimensions (UED). We will follow the well 
known procedure in literature. 11 

The relic abundance of dark matter \ is found by solving the Boltzmann 
equation for the evolution of the \ number density n 



where H is the Hubble parameter, v is the relative velocity between two 
x's, (o~v) is the thermally averaged total annihilation cross-section times 
relative velocity, and n eq is the equilibrium number density, (av) is often 
approximated by the non-relativistic expansion 



where x = ^ is the ratio of the dark matter mass to the temperature. By 
solving the Boltzmann equation analytically with appropriate approxima- 
tions, 11 the abundance of \ is given by 



^ = -3ffn-M(n 2 -<), 



(1) 



(av) = a + b(v 2 ) + 0((v 4 )) « a + 6b/x + O 




n x h 2 



1.04 x 10 9 x F 1 



(3) 



Mpi ^Jg*(x F ) a + 36/ x F 
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where the Planck mass Mpi = 1.22 x 10 19 GeV and <?* is the total number 
of effectively massless degrees of freedom. The freeze-out temperature, Xf, 
is found iteratively from 



where the constant c is determined empirically by comparing to numerical 
solutions of the Boltzmann equation. At the end the calculation of the relic 
abundance becomes computation of annihilation cross sections of relevant 
processes. 

In the case of Minimal UED (MUED), KK photon is a good dark matter 
candidate and a pair of KK photons annihilates to the SM fermions and 
Higgses. Here are the steps to get annihilation cross sections for CalcHEP . 

(1) Download models files for the dark matter model, in this case MUED 



(2) Unpack: tar -xvf MUED. tar and import this model by selecting "IM- 
PORT of MODELS" in the main CalcHEP menu (One can simply copy 
files into models under user's working directory but the model number 
needs to be changed consecutively). 

(3) Run Bl Bl — > e E for a typical KK photon annihilation. Users should 
see 4 diagrams with two kinds of KK fermions. 

(4) After squaring diagrams, go to "Symbolic calculations" and save the 
squared matrix elements in Mathematica code. Users should see that 
"symbl.m" is generated under "results" directory. 

(5) Now in Mathematica, run "sum_int.m" which is found in a directory 
called "utile" (or slightly modified version can be downloaded from 
http://susy.phsx.ku.edu/-kckong/tasi/TASIJVIUED_BlBl_ee.nb.). 

(6) After loading the package, run the following command in Mathematica. 

sum = 

addToSum:= sum 

= sum + Simplify [ totFactor numerator/denominator 
/. substitutions ]; 



e The model files can be obtained from the web site for this TASI tutorial, or directly 
from http://susy.phsx. ku.edu/~kckong/tasi/MUED. tar. 12 It is important to check the 
model file has no errors with CalcHEP version that users are running. This can be done 
by simply selecting "CHECK MODEL" under "Edit model" in the main menu. 





e 
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(7) Then load the generated matrix element in "symbl.m". By simply typ- 
ing "sum" we obtain the squared matrix element, \M\ 2 . "symbl.m" in- 
cludes squared matrix elements as well as the corresponding diagrams. 
Each squared diagram is computed by multiplying three quantities: 
totFactor, numerator and denominator. 

(8) Now the corresponding annihilation cross section is obtained by inte- 
grating over phase space. 



1 



16ttAi2 j to 



t-n 

2 



dt\M\ 2 , (5) 



where Ajj = X(s,m 2 ,m 2 ) = (s — m? — m 2 ) 2 — 4m 2 m 2 for a particle 
scattering, 1,2 — > 3,4. The lower (upper) bound is given by to(t w ) = 
((ml — m\— m% + m 2 ) 2 — (VA12 T VA34)) . Here s, t and u are the 
Mandelstam variables and mi = m = m 2 and m 3 = = m 4 in this 
example. 

(9) By simplifying the cross section, we should obtain the following expres- 
sion 

a(BlBl -> e+e") = ^j^ns 2 ^ ( 10 ( 2m2 + ~ 7 P S ) > 
'"" (6) 



where /3 = y 1 — F's are the Hypercharges of electrons and 51 is 
the strength of the U(1)y interaction. 
(10) Finally two leading terms (a and b terms) arc obtained by expanding 
in terms of the relative velocity. 



2.4. Implementing new particles and new interactions 
2.4.1. Higgs effective couplings 

For the Higgs effective couplings to photons and gluons, we do not need 
to include new particles but need to include new interactions of the Higgs 
(new in a sense that CalcHEP does not have it). For a relatively light Higgs, 
these may be expressed in terms of dimension 5 operators, 

tggh = -^9gghG c l v G^ a h 1 (7) 
-£77/1 — — ~^9^~yhFj_ii/F^ i (8) 
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and the couplings are calculated approximately as 

( 7 2 2 26 , \ 
' 9 ^-^( v 1+ 30 ++ 2i r * + 525 T ' + ---J' (9) 

_iilZ (^ — 2 3 5248 4 

9llh ~ ™ 18 ^ + 235^ + 1645 Ttu + 8225 ^ + 90475^ 

1280 5 54528 6 56 32 2 \ 

+ 29939 Tl " + 1646645 ^ ~ 705 T ' ~ 987 T * + " ' J ' ( ^ 

2 2 

where r f = and r w = ■ One can compnte Feynman rules for 

these interaction easily and implement them in "Vertices" in CalcHEP s . 

Standard Model 
Vertices 

Al |A2 I A3 |A4 |> Factor <|> Lorentz part <| 



G 


IG 


Ih 1 


1 (GG-2/(4*pi))/(12*pi*vev)*(-4) Ipl .p2*ml .m2 - 


m2.pl*ml.p2 


A 


IA 


Ih 1 


IchAA Ipl.p2*ml.m2 - 


m2.pl*ml.p2 



We only take the leading term in the couplings for simplicity. The struc- 
ture of this file is very easy to understand, which is again big advantage of 
using CalcHEP . In the "Lorentz" part, "pi" denotes the momentum of the 
first particle in the particle listing, Al, which is G, "p2" is the momentum 
of the second particle, A2, which is G again in this case, "ml" and "m2" 
are Lorentz indices and hence pl.m2 (or m2.pl) means the momentum of 
the first particle Al, which carries the Lorentz index of the second particle 
A2. The ml.m2 is the metric, g m im2- In the "Factor" we can include full 
expression of the coupling as we have done for g gg h or we can define the 
coupling in the "Constraints" , 

Standard Model 
Constraints 

>Name <|> Expression <| 



chAA |hAA(EE,vev) 



and include definition of hAA(EE,vev) in the "usrfun.c" . 



^MadGraph 13 has the same implementation. See the following wcbpagc for more informa- 
tion: https://server06.fynu.ucl.ac.be/projects/madgraph/wiki/Models/HiggsEffective. 
g For a different approach, users are encouraged to investigate FcynRulcs 14 and Lan- 
HEP. 15 
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double hAA (double EEE, double VW) 
{ 

double hAAcoupling; 
double alpha; 
double pi; 

pi=acos(-l) ; 

alpha = EEE*EEE/ (4*pi) ; 

hAAcoupling = alpha/ (pi*VVV) * 47/18; 
/* first term from the effective coupling of 

the higgs to two photons*/ 

return hAAcoupling; 

} 

One of the mode files "extlib#.mdl" has path to this user-defined code. 

SM with Gprime 
Libraries 

External libraries and citation <| 
$CALCHEP/utile/usrfun.c 



2.4.2. Implementing a color-octet vector boson 

A new particle can be implemented easily in CalcHEP using GUI interface. 
Suppose we are interested in a model with one new particle, color-octet 
vector boson G', which has interaction with the SM quarks (q and q), 

where A a is the Gell-Mann matrix, g% is the coupling strength of QCD and 
7 M s are the Dirac 7 matrices. 

(1) First we define the particle in "Particles". 

SM with Gprime 
Particles 

Full name I >A < I >A < I number I 2*spin I mass I width I color I aux I >LaTex(A) < I >LaTeX(A+) < I 

gluon IG IG 121 12 10 10 18 IG Ig Ig 



t-quark It IT 16 II IMt Iwt 13 I It I \bar{t} 

Gprime I ~G I ~G 1310002112 IMGP MwGP 18 IG I I 
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Most of necessary inputs are easy to understand. The exclamation 
mark "!" in the particle width means CalcHEP will run the decay pro- 
cess and calculate the width on the fly. The newly introduced variable, 
MGP is the mass of the new particle and we will define it in "Param- 
eter" . 

(2) For the interactions, CalcHEP already knows about the Gell-Mann 
matrices from the quantum charges and users do not need to include 
them in the vertices. 

SM with Gprime 
Vertices 



Al 


IA2 


IA3 


IA4 l> 


Factor 


< 1 > Lorentz part 




D 


Id 


|-G 


1 1 GG 




IG(m3)*(qV + qA * 


G5) 


U 


111 


|-G 


1 1 GG 




IG(m3)*(qV + qA * 


G5) 


T 


It 


|-G 


1 1 GG 




IG(m3)*(tV + tA * 


G5) 



The "G" and "G5" denote the Dirac 7 matrices and the "GG" is the 
strong coupling constant. Since we have introduced new variables (qV, 
qA, tV, tA and MGP), we should define them either in "Parameters" 
or in "Constraints" . We do this in the "Parameters" for this exercise. 

SM with Gprime 
Parameters 



>Name < I Value I > Comment < I 

alfEMZ 10.0078180608 IMS-BAR electromagnetic alpha(MZ) 



qA II I axial coupling of Gprime to qqbar 

qV 10 I vector coupling of Gprime to qqbar 

tA |-1 I axial coupling of Gprime to ttbar 

tV 10 I vector coupling of Gprime to ttbar 

MGP 1 1000 I mass of Gprime 



This completes implementation of a color-octet vector boson and its 
interaction to the SM quarks. 

3. PYTHIA 

PYTHIA is frequently used for event generation in high-energy physics. The 
emphasis is on multi-particle production in collisions between elementary 
particles. This in particular means hard interactions in e+e~, pp and ep 
colliders, although also other applications are envisaged. The program is 
intended to generate complete events, in as much detail as experimentally 
observable ones, within the bounds of our current understanding of the 
underlying physics. 

In this tutorial, we will consider a few specific examples. Unlike CalcHEP 
, PYTHIA requires a main driver, which will be compiled together with 
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PYTHIA code. We will use Fortran as our main compiler 11 . 
3.1. Installation 

Make sure you have a Fortran compiler, e.g. g77, gfortran and 
etc. We will use PYTHIA version 6.4.25 and it can be down- 
loaded from http: //home, thep . lu. se/^torbjorn/Pythia. html. Create 
"pythia-6. 4.25.o" by typing "g77 -c pythia-6.4.25.f" . You should see pythia- 
6.4.25.0 if the compile was successful. If you have a different compiler, 
use that one instead. Download an example Fortran code ("main61.f") for 
the most trivial test, and copy it to the same directory where you have 
PYTHIA . Compile it with PYTHIA and make an executable, by typing "g77 
-o main61.x main61.f pythia-6. 4.25. o" . Run "./main61.x" in your terminal 
and take a look at output on your screen. 

To make compile processes easier, often "makefile" is used. An example 
is shown below. 

# Makefile for pythia examples 

FF = gfortran 

# FF = g77 

# this is ok with gfortran from fink 

FFLAGS = -g -static -w -f no-second-underscore 

# For MAC: use the following FFLAGS 

# FFLAGS = -gdwarf-2 -static -w -f no-second-underscore 

DBJS = main61.o . ./pythia-6. 4.25.0 
EXEC = pythia. x 

#LIB = /cern/pro/lib/libmathlib . a 

all: $(EXEC) 

$(EXEC) : $(0BJS) 

$(FF) -o $(EXEC) $(0BJS) $(LIB) 

.f .0: 



h We will assume that users are familiar with Fortran. The most recent version is PYTHIA 
8 and is written in C++. 
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$(FF) $(FFLAGS) -c $< 

In this example of makefile, the executable is named as "pythia.x" and 
gfortran is used as a Fortran compiler. If additional libraries are needed, one 
can add them to the "LIB" variable and uncomment out the corresponding 
line. 

3.2. Running examples 
3.2.1. ti production 

The beginning of PYTHIA main driver include many lines of common blocks 
and we will not go over their details. Instead we focus on specific examples. 
The first example with PYTHIA is ti production at the Tevatron. Here we 
show general structure of the main code 1 . It is readable, if users understand 
basics of collider physics. 

c 

C. . .First section: initialization. 

c 

C. . .Number of events 
nevpythia=10000 

C. . .MSEL selects individual process 
MSEL=6 ! t quark, 

C... Sample code to force only decays of interest 
do i=190,208 ! Turn off all W decays 

mdme(i,l) = ! See the manual for the meaning of mdme 
enddo 

mdme (206,1) = 1 ! Turn on electron + neutrino 

C particles masses 

PMAS(6,1)=172. ! top mass 

PMAS(24,1)=80. ! W mass 

C...If interested only in cross sections and resonance decays: 
C... switch off initial and final state radiation, 



'The code below is not complete and users will not be able to compile with this. 
Part of the code is shown for discussion. Complete source can be found from 
http://susy.phsx.ku.edu/~kckong/tasi/. It also contains many comments. They do not 
affect the running of the code and are added for user's convenience. 
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C. . .multiple interactions and hadronization. 

MSTP(61)=0 ! initial state radiation [0] off, [D=2] on 

[1] : on for QCD radiation in hadronic events 

and QED radiation in leptonic ones 
[2] : on for QCD/QED radiation in hadronic events 
and QED radiation in leptonic ones 
MSTP(71)=0 ! final state radiation [0] off, [D=l] on 
MSTP(81)=0 ! multiple interactions [0] off, [D=l] on 
MSTP(111)=0! fragmentation and decay [0] off, [D=l] on 
MSTP(91)=0 ! No primordial kT 



C. . 

c 
c 
c 
c 



PDF: [D=8] , page 200 



MSTP(32) = 4 
MSTP(32) =11 
MSTP(52)=2 
MSTP(51)=7 



Q~2 value set, root-s-hat 

Q~2 = (m3+m4) ~2/4 

to use external PDF lib, [D=l] 

[D=7] 

choose pdf. 7:CTEQ5L L0, 8: CTEQ5M1 



C. .. Initialization for the Tevatron or LHC or ILC. 



CALL PYINIT ( ' CMS ' , ' e+ ' , ' e- ' , 500D0) 
CALL PYINIT ('CMS' , 'p' , 'pbar' ,1960D0) 
CALL PYINIT ( ' CMS ' , ' p ' , ' p ' , 14000D0) 
CALL PYINIT ( ' CMS ' , ' p ' , ' p ' , 7000D0) 



! 500 GeV ILC 
! Tevatron 
! 14 TeV LHC 
! 7 TeV LHC 



C initialize histograms 

CALL PYB00K(1 , 'ttbar invariant mass', 
CALL PYB00K(2,'PT distribution of top', 
CALL PYB00K(3, 'inv mass of b and e-', 
CALL PYB00K(4, 'inv mass of bbar and e-'. 
CALL PYB00K( 5, 'minimum sqrt(s)', 



100,0D0,1000D0) 
100,0D0,1000D0) 
50.0D0, 300D0) 
50.0D0, 300D0) 
100,0D0,1000D0) 



C... Second section: event loop. 

c 

C. . .Loop over the number of events. 

DO IEV=l,nevpythia 

IF(M0D(IEV,100) .EQ.0) print*, 'Now at event number', IEV 

C... Event generation. 
CALL PYEVNT 

CALL PYHEPC(l) ! convert to HEPEVT standard 

! see section 5.4 of the manual 



C...List first few events, say 2. 
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IFCIEV.LE.2) then 

CALL PYLIST(l) 

read* 
END IF 

c... Example histogram: look at the ttbar invariant mass distribution 

do ihep=l,nhep ! loop through all particles in the event record 
if (isthep(ihep) . eq. 3) then ! only look at the summary portion 



if ( 


idhep(ihep) 


.eq. 


6 ) 


it 




ihep 


! found 


top 


if ( 


idhep(ihep) 


.eq. 


-6 ) 


itbar 




ihep 


! found 


anti-top 


if ( 


idhep(ihep) 


.eq. 


-11) 


ie 




ihep 


! f ound 


e- 


if ( 


idhep(ihep) 


.eq. 


+ 11) 


iebar 




ihep 


! found 


e+ 


if ( 


idhep(ihep) 


.eq. 


5 ) 


ib 




ihep 


! found 


b 


if( 


idhep(ihep) 


.eq. 


-5 ) 


ibbar 




ihep 


! found 


bbar 


if( 


idhep(ihep) 


.eq. 


12) 


in 




ihep 


! found 


nu_e 


if ( 


idhep(ihep) 


.eq. 


-12) 


inbar 




ihep 


! found 


nu_ebar 



endif 

enddo ! loop through all particles in the event record 



C Fill histograms and calculate things. 

CALL PYFILL(1, inv_mass (it , itbar) , 1D0) 
CALL PYFILL(2, pt(it), 1D0) 
CALL PYFILL(3, inv_mass (ie , ib) , 1D0) 
CALL PYFILL(4, inv_mass (ie , ibbar) , 1D0) 

cc Smin 

do i=l,4 

pvis(i) = phep(i , ie)+phep(i , iebar)+phep(i , ib)+phep(i , ibbar) 
enddo ! the total momentum of all visible particles 

ptvis = dsqrt( pvis(l)**2d0 + pvis(2)**2d0 ) 

! PT sum of all visible particles 
mvis = pvis(4)**2d0 - pvis(l)**2d0 - pvis(2)**2d0 - pvis(3)**2d0 

! visible mass 
if (mvis .It. OdO) mvis=0d0 
mvis = dsqrt(mvis) 

smin = dsqrt( ptvis**2d0 + mvis**2d0 ) + ptvis 

CALL PYFILL(5, smin, 1D0) 

c Event selection of a CDF analysis for A_FB in dilepton channel 
c 

c CDF Note 10436, 5.1 fb-1 

c A_fb =0.42 0.15stat 0.05syst 



! wrong pair 
! correct pair 
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c A_fb(theory) =0.06 0.01 
c 

c electrons : CAL ET > 20 GeV (little energy in HCAL) , 
c I eta I < 1.1 and 1.2 < I eta I < 2.8 

c muons: track PT > 20 GeV, I eta I < 1.0 
c MET > 25 GeV 

c jets: PT > 25 GeV, I eta I < 2.5 

c HT > 200 GeV (HT scalar sum of MET, leptons, jets) 

MPT(l) = phep(l,in)+phep(l,inbar) 
MPT (2) = phep(2,in)+phep(2,inbar) 
MET = dsqrt( MPT(l)**2dO + MPT(2)**2dO ) 



HT = MET + pt(ie) + pt(iebar) + pt(ib) + pt(ibbar) ! ignoring ISR 

if( abs(eta(ie)) .It. 2. 8 .and. pt(ie) .gt.20d0 ! e- 

& .and. abs(eta(iebar)) .lt.2.8 .and. pt (iebar) .gt . 20d0 ! e+ 

& .and. abs(eta(ib)) .It. 2. 5 .and. pt(ib) .gt.25d0 ! b 

& .and. abs(eta(ibbar)) .lt.2.5 .and. pt (ibbar) . gt . 25d0 ! bbar 

& .and. MET .gt.25d0 ! MET 

& .and. HT .gt.200d0 ! HT 



& ) then 

counter = counter + 1 

end if 

ENDD0 ! Loop over the number of events. 

c 

C. . .Third section: produce output and end. 
c 

C... Cross section table. 
CALL PYSTAT(l) 

C. . . Finalize analysis and report results 
C. . . Plot histograms. 

c PYDUMP(MDUMP,LFN,NHI,IHI) 
c MDUMP=3: (x,y) format 
c LFN: file number 

c NHI: number of histograms to be dumped; 
c if then all existing histograms are dumped, 

c IHI : array containing histogram numbers 
c in the first NHI positions for NHI nonzero 

CALL PYDUMP(3, 1,1,1) ! ttbar inv mass 
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CALL PYDUMP(3,2,1,2) 
CALL PYDUMP(3,3,1,3) 
CALL PYDUMP(3,4,1,4) 
CALL PYDUMP(3,7,1,5) 



PT of the top 

(e-,b) - wrong pair 

(e-,bbar) - correct pair 

smin 



print*, 'Total number of events before cuts 
print*, 'Total number of events after cuts 
print*, 'efficiency 



nevpythia 
counter 



& 



dble (counter ) / dble (nevpythia) 



The main code is divided into 3 sections; initialization, event loop and 
ending. In the first part, users set up process, the number of events, collider 
type, center of mass energy, beam environments, masses, decay patterns, 
initializing histograms etc. The second part is for the event generation, 
where most of actual analysis is done. Users are supposed to store relevant 
information of each event, since PYTHIA does not save event information 
by default. Users fill histograms in this second part. We will use PYTHIA 
commands for plots but often many users use their own favorite plotting 
program and do not use plotting package in PYTHIA . In the last section, 
histograms may be plotted, users can finalize analysis and report results. 
We will go over a few important commands below. 

(1) First of all, we save the number of events to generate in a variable, 
nevpythia. Then we use MSEL variable to choose a process to be sim- 
ulated. Selected other processes are shown below. 



.Standard Model 
MSUB(l)=l 
MSTP(43)=2 



MSUB(18)=1 

CKIN(l)=100dO 

MSUB(22)=1 

MSUB(23)=1 

MSUB(24)=1 

MSUB(25)=1 

MSUB(26)=1 

MSUB(27)=1 



Drell-Yan 

[1] only photon diagram included 
[2] only Z diagram included 
[3] both diagrams + interference 
f fbar -> gamma gamma 

f fbar -> Z Z 

f fbar -> W Z 

f fbar -> Z h 

f fbar -> W+ W- 

f fbar -> W h 

f fbar -> h h 



C. 

c 



.New physics: Doubly charged Higgs 

MSUB(349)=1 ! f fbar -> HL++ HL— 
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c MSUB(350)=1 ! f fbar -> HR++ HR— 

c MSUB(351)=1 Iff -> f f HL 

c MSUB(352)=1 Iff -> f f HR 

c kc = pycomp(9900041) 

c PMAS(kc,l)=1000 

c kc = pycomp (9900042) 

c PMAS(kc,l)=1000 

c. . .top quark 

MSEL=6 ! t quark 

C. . .New physics: fourth generation/Little Higgs 
c. . .t> 

c MSEL=8 ! fourth generation t' pair production 

c MSTP(7)=8 ! choice of heavy flavor, superseeded by MSEL=4-8 

c PMAS(8,1)=500 ! mass of the t> 

c MSEL=38 ! fourth generation single t' 

c MSUB(83)=1 ! q f -> Q f 

c. . .b' 

c MSEL=7 ! fourth generation b' pair production 

c MSTP(7)=7 ! choice of heavy flavor, superseeded by MSEL=4-8 

c PMAS(7, 1)=400 ! mass of the b> 

c MSEL=37 ! fourth generation single b' 

c MSUB(83)=1 ! q f -> Q f 

c MSTP(1)=4 ! four generations 

c do i=56,75 ! Turn on b' and t' decays 

c mdme(i,l) = 1 

c enddo 

c MSTP(127)=1 ! do not crash if vanishing cross-sections 



(2) Users are allowed to force a particular decay of certain particle, W 
in this example. This is done with "mdme" switch. To list all decay 
channels, widths etc, one can use "CALL PYSTAT(2)". 

(3) PMAS sets masses. 

(4) Initial and final state radiation, multiple interactions and fragmcnta- 
tion/hadronization are treated by MSTP switches. 

(5) To initialize the collider type, use "PYINIT", for instance, "CALL 
PYINIT(CMS,p,pbar,1960D0)" for Tevatron. 

(6) Event generation is performed by calling "PYEVNT" and information 
can saved into HEPEVT standard format for further analysis. 

(7) "CALL PYLIST(l)" shows basic event listing, which is self- 
explanatory. 

Event listing (summary) 
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I particle/jet KS KF orig p_x p_y p_z E m 



1 


!p+! 


21 


2212 








000 





000 


980 


000 


980 


000 





938 


2 


!pbar-! 


21 


-2212 








000 





000 


-980 


000 


980 


000 





938 


3 


!u! 


21 


2 


1 





000 





000 


479 


080 


479 


080 





000 


4 


! ubar ! 


21 


-2 


2 





000 





000 


-137 


794 


137 


794 





000 


5 


!u! 


21 


2 


3 





000 





000 


479 


080 


479 


080 





000 


6 


! ubar ! 


21 


-2 


4 





000 





000 


-137 


794 


137 


794 





000 


7 


!t! 


21 


6 





-39 


462 


145 


528 


32 


386 


232 


316 


173 


741 


8 


! tbar ! 


21 


-6 





39 


462 


-145 


528 


308 


899 


384 


558 


172 


426 


9 


!W+! 


21 


24 


7 


-65 


111 


123 


272 


76 


449 


176 


479 


76 


583 


10 


!b! 


21 


5 


7 


25 


648 


22 


256 


-44 


063 


55 


837 


4 


800 


11 


!W-! 


21 


-24 


8 


77 


947 


-79 


254 


123 


668 


184 


042 


78 


871 


12 


! bbar ! 


2 


-5 


8 


-38 


485 


-66 


274 


185 


231 


200 


517 


4 


800 


13 


!e+! 


21 


-11 


9 


-3 


487 


85 


329 


30 


144 


90 


564 





001 


14 


! nu_e ! 


21 


12 


9 


-61 


624 


37 


943 


46 


305 


85 


915 





000 


15 


!e-! 


21 


11 


11 


69 


103 


-30 


988 


109 


185 


132 


879 





001 


16 


!nu_ebar ! 


21 


-12 


11 


8 


844 


-48 


266 


14 


483 


51 


163 





000 



(8) It is highly recommend to save particle IDs in some variables as shown 
in the example. This is done using "idhep" command and PDG num- 
bers. 

(9) Users can compute certain quantities and do analysis afterwards. In 
this example, we are making 5 histograms, invariant mass of the top 
pair, transverse momentum of the top, invariant mass of a lepton and 
b/b, and y/s min . 

(10) In the last section of the main code, users can produce output/figures 
and end the program. For instance, the output of this example is shown 
below. The cross section table is called by using "CALL PYSTAT(l)". 
The second part is efficiency due to some cuts that are introduced. It 
shows that 51% of total 10K events passed those cuts. 

********* PYSTAT : Statistics on Number of Events and Cross-sections ********* 



I I II 

I Subprocess I Number of points I Sigma I 

I I II 

j 1 1 (mb) I 

I I II 

I N : o Type I Generated Tried I I 

I I II 



I I II 

I All included subprocesses I 10000 201243 I 7.174D-11 I 

I 81 q + qbar -> Q + Qbar, mass I 9402 185523 I 6.757D-11 I 

I 82 g + g -> Q + qbar, massive I 598 15720 I 4.170D-12 I 

I I II 



********* Total number of errors, excluding junctions = ************* 

********* Total number of errors, including junctions = ************* 

********* Total number of warnings = ************* 

********* Fraction of events that fail fragmentation cuts = 0.00000 ********* 



August 2, 2012 0:9 



WSPC - Proceedings Trim Size: 9in x 6in CalcHEP-PYTHIA 




Total number of events before cuts: 10000 
Total number of events after cuts : 5114 
efficiency : 0.51139999999999997 



The other outputs are figures. This example uses PYTHIA commands 
("PYDUMP"), to write results into files. 

With help of a plotting program, we show the following results. Fig. 2 
shows various distributions in the tt dilepton production: tt invariant mass 
(a) , transverse momentum of the top (b) , invariant mass of a b-quark and 
a lepton (c), and \fs min (d). In Fig. 2(c), both wrong and correct combi- 
nations are shown, which are different and therefore can be used to reduce 
combinatorial background. 16 Even in the case of two missing neutrinos the 
y/s min exhibits a peak at 2Mj, which is shown in the vertical line. 17,18 
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3.2.2. slepton production 

In this example, we consider pair production of the right handed selectron 
at an 500 GeV ILC. 

(1) Main structure of the code is very similar to the previous exam- 
ple. The process is again set up by "MSEL". In this case we use 
"MSUB(202)=1". Other examples of SUSY processes are shown be- 
low and a complete list is found in the manual. Users have various 
options to dehnc mass spectrum for SUSY production: PYTHIA RG 
running, external mass spectrum generators or LHA input files. We 
will consider one of mSugra point, SPSla, and use PYTHIA commands, 
"IMSS" and "RMSS", which define mass spectrum. 



.Generic SUSY simulation (SUGRA scenario) 



MSEL=39 
MSEL=42 
MSEL=0 
MSUB(202)=1 
MSUB(201)=1 
MSUB(204)=1 
MSUB(208)=1 

MSUB(210)=1 
MSUB(212)=1 

MSUB(213)=1 
MSUB(214)=1 

IMSS (1) =2 

IMSS(1)=12 
RMSS(1)=250D0 
RMSS(4)=1D0 
RMSS(5)=10D0 
RMSS(8)=100D0 
RMSS(16)=0D0 



! All MSSM processes at once 

! Slepton production 

! One by one 
! f fbar -> ~e_R ~e_Rbar 

! f fbar -> ~e\_L ~e\_Lbar 

! f fbar -> ~mu\_L ~mu\_Lbar 

! f fbar -> ~tau2 ~tau2bar 

! f fbar -> ~1\_L ~nu\_Lbar 
! f fbar -> ~tau\_2 "nutaubar 

! f fbar -> "mil "nulbar 

! f fbar -> "nutau "nutaubar 

! MSUGRA spectra from analytic approximation 
! MSUGRA spectra from IsaJet 
Mhalf : common gaugino mass 

Sign of mu (magnitude irrelevant if IMSS(1)=2) 
tan beta 

MO: common scalar mass 
AO: common trilinear coupling 



(2) The collider type and energy are set by "PYINIT" . 



CALL PYINIT ( ' CMS ' , ' e+ ' , ' e- ' , 500D0) 



! 500 GeV ILC 



(3) The rest is basically the same as all other examples. Initializing his- 
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tograms, identifying particles, calculating relevant quantities, filling 
histograms, calculating efficiency, plotting histograms, report results 
and etc. For instance, to calculate Mt2, one needs momentum of visi- 
ble particles and the missing transverse momentum. In this case, two 
visible particles are electron and positron, and the missing momentum 
is balanced by the momentum of e + e~ pair. 

do ihep=l,nhep ! loop through all particles in the event record 
if (isthep(ihep) . eq. 3) then ! only look at the summary portion 
if( idhep(ihep) .eq. -11) ie = ihep ! found e- 
if( idhep(ihep) .eq. 11) iebar = ihep ! found e+ 

endif 

enddo ! loop through all particles in the event record 

pinl(l) = phep(l.ie) 
pinl(2) = phep(2,ie) 
pin2(l) = phep(l , iebar) 
pin2(2) = phep(2,iebar) 

mchi = PMAS (310,1) ! KC = PYCDMP (1000022) 
mt2 = mt2_sleptons(pinl,pin2,mchi) 

The four momenta of particles are saved in "PHEP" variable. 

Various kinematic distributions are shown in Fig. 3: the invariant mass 
of the e+e~ (a), Mt2 (b), the energy distribution (c) and the transverse 
momentum (d) of the electron. The energy distribution and the Mt2 exhibit 
interesting end point structures which are useful for extracting particles 
masses, sleptons and neutralinos in this example. 19 ' 20 



3.2.3. p T -^ttt + W 

In this example, we simulate techni p production with PYTHIA at the Teva- 
tron, setting masses M p o = 290 GeV and M^± — 160 GeV. We consider 

the p\ decay to iri^ + W T and 7r^ to 2 jets. For faster event generation, we 
force the decay of W into electron and neutrino. 

C. . .Possibility to set masses freely: 
C. . .pi_tech0 

PMAS(PYC0MP(KTECHN+111) ,1)=160D0 ! KTECHN=3000000 
C. . .pi_tech+-. ! for techi-particles 

PMAS(PYC0MP(KTECHN+211) ,1)=160D0 
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50 100 150 200 350 300 

M e .„- (GeV) 



120 140 

M T2 (GeV) 




Fig. 3. Various distributions for selectron production at a 500 GeV ILC. (a) M e + £ 
(b) M T2 (c) E e (d) P«. 



C. . .rho_tech0 

PMAS (PYC0MP (KTECHN+1 13) , 1) =290D0 

C. . .rho_tech0. 

MSUB(191)=1 
C . . . rho_tech+- . 
c MSUB(192)=1 
C. . .omega_tech. 
c MSUB(193)=1 



.Sample code to force only decays of interest 
do 1=190,208 

mdme(i,l) = ! Turn off all W decays 
enddo 
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100 200 300 400 100 200 300 400 

m wjj (GeV) sqrt(s ml „) (GeV) 

Fig. 4. Invariant mass distribution (a) and y/s min (b) for pJJ, — > ir^ + W^ — > jj + Pfe 
at the Tevatron. 



mdme (206 , 1) = 1 ! Turn on electron + neutrino 

Forming an invariant mass of the two jets, one can easily find a bump. 
With fully reconstructed W, the p% appears as a resonance in the invariant 
mass of W + jj, which is shown as the blue histogram in Fig. 4(a). In 
the leptonic decay of the W, the transverse momentum of the neutrino 
is determined by the missing transverse momentum, assuming that the 
neutrino is the only missing particle in the event. Mass-shell condition of 
the W provides the z-componcnt of the neutrino momentum up to two fold 
ambiguity. Invariant mass of two jets and tvn system with both solutions 
is shown in red. With appropriate detector effects, the invariant mass get 
further smeared as shown in back in Fig. 4(a). Without reconstructing the 
neutrino momentum, one can obtain mass information looking at the end 
point structure in the \fs min distribution as shown in Fig. 4(b). The red 
vertical lines shows location of p9p resonance. 

3.2.4. The same sign dilepton in SUSY: LM6 

We consider a pair production of gluino for LM6 point (a mSugra point in 
CMS) at the 14 TeV LHC. This time we define mass spectrum from a file in 
LHA format, which is generated by SuSpect. 21 This is done with "IMSS" 
and "PYSLHA" , as shown below. 

MSUB(244) =1 ! g g -> gluino gluino 

c 

c IMSS(1)=2 ! MSUGRA spectra from analytic approximation 

cc IMSS(1)=12 ! MSUGRA spectra from IsaJet 

c RMSS(l) =400D0 ! Mhalf : common gaugino mass 
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c RMSS(4) =1D0 ! Sign of mu (magnitude irrelevant if IMSS(1)=2) 

c RMSS(5) =10D0 ! tan beta 

c RMSS(8) =85D0 ! MO: common scalar mass 

c RMSS(16)=0D0 ! AO: common trilinear coupling 



C. 
C. 
C. 

c. 
c. 
c. 
c. 
c. 



open(99 , FILE= ' suspect . txt ' ) 
IMSS(l) =11 ! for SLHA input 

IMSSC13) = ! 0=MSSM particle content, 1=NMSSM 

IMSS(21) = 99 ! Logical Unit Number for SLHA spectrum read-in. 

IMSS(22) = 99 ! Logical Unit Number for SLHA decay read-in. 

! normally this should go together for consistency 
CALL PYSLHA(5,0,IFAIL) 
CALL PYSLHA(2,0,IFAIL) 
close(99) 



MUPDA=0 
MUPDA=1 
MUPDA=2 

MUPDA=3 
MUPDA=4 
MUPDA=5 



READ Q NUMBERS/PARTICLE ON LUN=IMSS(21) 

READ SLHA SPECTRUM ON LUN=IMSS(21) 

LOOK FOR DECAY TABLE FOR KF=KF0RIG ON LUN=IMSS(22) 

(KF0RIG=0 : read all decay tables) 

WRITE SPECTRUM ON LUN=IMSS(23) 

WRITE DECAY TABLE FOR KF=KF0RIG ON LUN=IMSS(24) 
READ MASS FOR KF=KF0RIG ONLY 
(KF0RIG=0 : read all MASS entries) 



c 


IMSS(l) 


= 1 ! A general 


MSSM simulation 


c 


RMSS( 1) 


= 100.0 


! M_l 500.0' 


c 


RMSS( 2) 


= 300.0 


M_2 


c 


RMSS( 3) 


= 1225.0 


M_3 


c 


RMSS( 4) 


= 1100.0 


! Mu 200.0' 


c 


RMSS( 5) 


= 10.0 


! Tan_beta 


c 


RMSS( 6) 


= 1800.0 


! M_S1_L 


c 


RMSS( 7) 


= 1800.0 


! M_S1_R 


c 


RMSS( 8) 


=5500.0 


M_Sq_L 


c 


RMSS( 9) 


=5500.0 


M_Sq_R 


c 


RMSS(IO) 


=5500.0 


M_Sq3_L 


c 


RMSS(ll) 


=5500.0 


M_Sbottom_R 


c 


RMSSC12) 


=5500.0 


M_Stop_R 


c 


RMSS(13) 


=5500.0 


M_Stau_L 


c 


RMSS(14) 


=5500.0 


M_Stau_R 


c 


RMSS(15) 


= 800.0 


A_b=Bottom trilinear coupling 


c 


RMSS(16) 


= 800.0 


A_t=Top trilinear coupling 


c 


RMSSC17) 


= 0.0 


A_tau=Tau trilinear coupling 


c 

c 


RMSS(19) 


= 400.0 


M_A=Psc higgs param. 



CALL PYINIT( 'CMS ' , 'p' , 'p' , 14000D0) ! 14 TeV LHC 
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Fig. 5. Mt2 (a) and invariant mass (b) distributions of the same-sign dilepton for LM6 
at the 14 TeV LHC. 
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Fig. 6. Missing transverse momentum distribution $ T for LM6. 



Fig. 5 shows the ID decomposed Mn (a) and invariant mass (b) dis- 
tributions of the same-sign dilepton for LM6 at the 14 TeV LHC and Fig. 
6 shows the distribution of the missing transverse momentum with and 
without detector effects. 



3.3. CalcHEP-PYTHIA Interface 

CalcHEP is often linked to PYTHIA for further simulation, especially to in- 
clude ISR, FSR, multiple parton interaction and etc. To use this feature, 
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we need an event-file in LHA format, which is generated by CalcHEP , and 
will be fed into PYTHIA . Suppose the file name is "event _mixer.lhe" , which 
we have generated in an earlier example, using "eventjnixer" . We also 
need "event2pyth.c" , which is in the CalcHEP "utile/" directory. 3 Here we 
should make sure a routine "UPEVNT" is deleted in this "event2pyth.c" , 
since PYTHIA already has one. The rest is to write a main code for PYTHIA 
, part of which is shown below. 

C...If interested only in cross sections and resonance decays: 
C... switch off initial and final state radiation, 
C. . .multiple interactions and hadronization. 



MSTP(61)=0 



MSTP(71)=0 
MSTP(81)=0 
MSTP(111)=1 
MSTP(91)=0 



initial state radiation [0] off , [D=2] on 

[1] : on for QCD radiation in hadronic events and 

QED radiation in leptonic ones 
[2] : on for QCD/Q.ED radiation in hadronic events and 

Q.ED radiation in leptonic ones 
final state radiation [0] off, [D=l] on 
multiple interactions [0] off, [D=l] on 
fragmentation and decay [0] off , [D=l] on 
No primordial kT 



c. . .USER MODE 

MSTP(161)=21 
MSTP(162)=21 

0PEN(2l, FILE= ' event_mixer . lhe ' , STATUS= ' UNKNOWN ' ) 
CALL PYINITC USER' , ' ',' ' ,0d0) 

Two Fortran files and one C file need to be compiled together. Make 
sure to include "event2pyth.o" in your makefile. 



DBJS = example_interface.o . . /pythia-6 . 4 . 25 . o event2pyth.o 



4. Summary 

In high energy physics there are various tools for different purposes. None 
of these tools is perfect and users should know their advantages and disad- 
vantages before choosing one for their study. We have looked at CalcHEP 
and PYTHIA among many tools. Both are commonly used event genera- 
tors. Especially CalcHEP is linked to micrOmegas for dark matter study 
and PYTHIA is a hardcore software for collider physics. Although we only 
caught a glimpse of what they can do, it is our hope that beginners get 
basic ideas behind complicated structures and are not afraid of using them, 
and more advanced users get usefulness out of specific examples. CalcHEP 
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and PYTHIA programs are continuously being improved and many current 
issues will be addressed in the future update. 
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